function fig2_17

alpha=1;
lambda=1;
%b=sqrt(2)*12;
b=16;
t=40;
nx=500;
x=linspace(0,230,nx);
for ix=1:nx
    k0=sqrt(x(ix)/(3*alpha*t));
    if x(ix)<40
        ix0=ix;
    elseif x(ix)<200
        ix1=ix;
    end
    asy(ix)=sqrt(2)*F(k0)*cos(k0*x(ix)-alpha*k0^3*t-pi/4)/sqrt(3*alpha*k0*t);
    h(ix)=sqrt(2/pi)*rhom2(x(ix),t,lambda,b);
end

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.2,'MR',0.01,'ML',0.065,'P',0.04)

hold on
plot(x,h,'b','LineWidth',1.5)
axis([0 210 -1 1])
plot(x(ix0:ix1),asy(ix0:ix1),'--r','LineWidth',1.5)
box on
grid on

set(gca,'ytick',[-1 -0.5 0 0.5 1])
set(gca,'yTickLabel',{'$-1$';'$-0.5$';'0';'0.5';'1'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

ylabel('$h(x,40)$','Interpreter','latex','FontSize',14)
xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)

legend(' Num',' Asy','Location','NorthEast','FontSize',13,'Interpreter','latex','AutoUpdate','off')

%exportgraphics(gcf,'/Users/mark/Desktop/kdv3.pdf')

function RR=rhom2(x,t,lambda,b)
%  Integrate  over  [a, b]  using composite simpson and romberg
%  S(k) = composite simpson;  R(k) = romberg with composite simpson
%  err = relative iterative error using romberg

a=0; b=lambda+5/b;

%  interv = number of subinternals
tol=1e-8;
n=8;
S(1)=simp(a,b,n,x,t);
err=10*tol;
R(1)=1;  errR(1)=-1;
k=1;
while err>tol
    k=k+1;
    n=2*n;
    S(k)=simp(a,b,n,x,t);
    R(k)=(16*S(k)-S(k-1))/15;
    err=abs(1-R(k)/R(k-1));
end
RR=R(k);


function y=F(k)
lambda=1; alpha=1; b=16;
y=(((exp(-(k-lambda)^2*b^2/4)+exp(-(k+lambda)^2*b^2/4))*b*sqrt(2))/4);

function y=f(k,x,t)
lambda=1; alpha=1; b=16;
y=F(k)*cos(-alpha*k^3*t+k*x);


%%% composite simpson: n = number of subintervals
function s=simp(a,b,n,x,t)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
s=f(xd(1),x,t);
for j=2:2:n
    ff2=f(xd(j+1),x,t);
    s=s+4*f(xd(j),x,t)+2*ff2;
end
s=h*(s-ff2)/3;














